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Abstract 

The cage model for polymer reptation is extended to simulate DC elec- 
trophoresis. The drift velocity u of a polymer with length L in an electric 
field with strength E shows three different regions: if the strength of field 
is small, the drift velocity scales as t" ~ E/L; for slightly larger strengths, 
it scales as v ~ i?^, independent of length; for high fields, but still E <^ 1, 
the drift velocity decreases exponentially to zero. The behaviour of the first 
two regions are in agreement with earlier reports on simulations of the Duke- 
Rubinstein model and with experimental work on DNA polymers in agarose 
gel. 
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I. INTRODUCTION 

A widely used tool to separate mixtures of DNA molecules by length is DC electrophore- 
sis. The DNA is confined to an agarose gel, and an electric field is applied. Since DNA is 
negatively charged, it moves towards the positive electrode as a result of this electric field. 
As the drift velocity depends on the length, DNA fragments with different lengths end up 
in different bands, and can therefore easily be separated. 

Since DNA fragments are usually much longer than the typical spacing between the gel 
strands, they are unable to move sideways. De Gennes ||I[] described the motion of a polymer 
in such an environment, and termed it reptation: the polymers move by diffusion of 'defects' 
along the chain of monomers. Each defect contracts the polymer by a certain amount of 
length, called its stored length. When a defect passes a monomer, the monomer is moved by 
this distance. Figure |I| shows an example where a defect travels past monomer B. 

Two models are widely used to simulate reptation: the repton model, introduced by 
Rubinstein [0], and the cage model, introduced by Evans and Edwards 0. In both models, 
monomers reside on sites of a simple cubic lattice (or in two dimensions a square lattice), 
and are connected by bonds; the dynamics consist of single-monomer moves. 

In the repton model, stored length consists of zero-length bonds. For this model, it was 
already proposed by Rubinstein 0, and later proven by Prahofer and Spohn |^ that the 
diffusion constant D of the polymer in the limit of long polymer length L obeys the scaling 
L^D = 1/3. For finite lengths, the diffusion constant is known numerically exact up to 
length 20 and from Monte Carlo simulations up to length 250 P|. The repton model has 
been adapted for the study of DC electrophoresis by Duke 0. This Rubinstein-Duke model 
has been studied numerically for lengths up to 400 . Simulations of this model are easy 
because it can, without loss of generality, be reduced to a one-dimensional model. 

In the cage model, stored length consists of a pair of anti-parallel nearest-neighbour 
bonds, called 'kinks'. The polymer diffusion constant in this model has been determined 
numerically exact for small L [H] , and with Monte Carlo simulations for polymers up to length 



200 P,p!0[|. As in the repton model, the polymer diffusion constant scales as D ~ L~^. In 
this paper, we extend the cage model to simulate a charged polymer in an electric field. 
Interestingly, we find differences in the scaling of the polymer drift velocity as compared to 
the Duke-Rubinstein model. 

In section |T| we describe the cage model and present how the model can be extended 



to simulate reptation in a non-zero electric field. We discuss in sections [II C| to pE| how 
efficient simulations can be achieved with multispin coding; these sections are not needed for 
understanding other parts of the paper. In section |T| we discuss the simulation approach 
in detail. In section |^ we discuss scaling arguments for the drift velocity. The results are 
presented in section which includes statements about the polymer shapes and comparisons 
to previous reports. 

II. CAGE MODEL 

The cage model describes a polymer of L monomers, located on the sites of an infinite 
cubic lattice. The monomers are connected by L — 1 bonds with a length of one lattice 
spacing. A single step of the Monte Carlo simulation consists of selecting randomly a 
monomer and, if it is free to move, moving it to a randomly selected location (possibly the 
current location). 

The monomers at both ends of the polymer are always free to move, but monomers in the 
interior of the polymer are only free to move when the two neighbours along the chain are 
located on the same adjacent lattice site. Other movements might result in an acceptable 
polymer configuration, but are ruled out because they would allow the polymer to move 
sideways, which is not reptation. One possible interior move is shown in figure 0. Every 
possible move occurs statistically with unit rate, setting the time scale. A single elementary 
move thus corresponds to a time increment of At = l/(2(iL), where d is the dimensionality 
of the lattice; in our case, d = 3. 



A. Electric field 

In solution, DNA becomes negatively charged with a fixed charge per unit length. We 
incorporate this into the cage model by assigning a negative charge q per monomer. The 
polymer is located in a homogeneous electric field E, that acts on these charged monomers. 

For two monomer positions ri and r*2; separated by a displacement ri2 = "^2 — ri, the 
difference in potential energy is given by f/ = qE ■ fi2. The ratio of the corresponding 
Boltzmann probabilities is 

where E = \E\ is the electric field strength, and r = E ■ fi2/\E\ is the displacement parallel 
to the field. 

In a Monte Carlo simulation, this ratio determines at which rates the monomers are to 
be moved along the field or against it. We choose the direction of the electric field along 
one of the body diagonals of the unit cubes, because then the x, y and z directions are 
equivalent, and within one elementary move, the displacement r takes only the two values 
±2/v^ times the lattice spacing. For convenience, the units are chosen in such a way that 
qr/ksT = ±1. 

Each monomer moves with a rate R^ = exp{E) for moves which lower the energy and 
R^ = exp{~E) for moves which raise the energy. When a monomer is selected and is able 
to move to lower (higher) energy, it will do so with a probability P^ (P), given by 

P+ = P~ = (2) 

d e^ + e-^' d e^ + e-^' ^ ' 

The time increment corresponding to one elementary Monte Carlo move is thus equal to 



B. Bond Representation 

As described in section the monomers are connected by bonds, where each bond has one 
of 2d possible orientations. One way of describing the polymer configuration is by specifying 
the location of the first monomer and the orientation of all bonds. The advantage of this 
notation is that only the position of one monomer has to be stored plus the orientations of 
all bonds. The polymer in figure ^, for example, is described by the position of the first 
monomer, on the left side of the figure, and +x-y+z+x-x+z. 

The dynamics can be described in terms of bonds. The bonds that are located on both 
ends of the polymer are always free to move. The internal bonds are free to move only when 
they are part of a pair of oppositely oriented neighbouring bonds (a kink). The first and 
last monomer in figure ^ can change to any new bond: +x, +y, +z, -x, -y or -z. The kink 
configuration +x-x can change in any new kink: +x-x, +y-y, +z-z, -x+x, -y+y or -z+z. 

C. Multispin Coding 

With multispin coding, many polymers can be simulated in parallel. We used an ap- 
proach similar to the one by Barkema and Krenzlin [^. The simulation that we performed 
used 64-bit unsigned integers to simulate 64 different polymers in parallel. As described in 
section [11 B , there are six directions a bond can point to, so each bond can be encoded using 



three binary digits. It is now possible to encode 64 bonds in three integers x, y and z, as 
shown in table |l[ 

In each iteration of the inner loop of the algorithm a random monomer i, < i < L, 
is selected. When an inner monomer is selected the two surrounding bonds are compared; 
if they are opposites, they are replaced by a randomly generated pair of opposite bonds. 
Section [II E| describes how to generate those bonds. The first and last monomers are special 
cases, which are described in section [II D| . 

To find the kinks in all of the 64 polymers, we use Equation (HI): 



A {y^-l®yi) (4) 

A {zi-i ®Zi). 
Monomer i is surrounded by bonds i — 1 and i. Bit j of ki is 1 if the surrounding bonds of 
monomer i of polymer j are in opposite directions. 

If a monomer can be moved, it will be relocated using a list of random kinks encoded in 
X, y and z. Bonds i — 1 and i that surround monomer i are replaced by x, y and z and their 
binary complements, respectively. Equation (|^) shows how this can be done: 

Xi-i = {^ki A Xj-i) V {ki A x) 

?/i_i = {^ki A yi_i) V (/cj A y) 

Zi-i = {^ki A Zi-i) y {kiAz) 

(5) 
Xj = {-^ki A Xj) V {ki A -ix) 

l/i = {-^ki A i/i) V (A;i A -^y) 

Zi = (-ifcj A Zj) V (fcj A -iz). 

These 27 logical operations replace the kinks near monomer i in all 64 polymers; polymers 

that have no kink near monomer i are left unaltered. 



D. First and last monomer 

The first and last monomers are always free to move. When one of those monomers is 
selected we can just replace the bonds with randomly generated bonds. To keep track of 
the position of the first monomer we have to calculate the distances traveled in the x, y and 
z directions. Since those directions are symmetric we only calculate r = x + y + z. For this 
we only need to know whether the first bonds point at a negative direction, which is one of 
-X, -y and -z. This is done using the following equation: 



d= {xqA yo) V {yo A Zq) V {zq A Xq) 



(6) 



Now the new random bonds are inserted: 



Xo = ->x 

yo = ^y C^) 

Zq = -^Z. 



Monomer only has one bond, which is number 0. Section |IIQ tells us that we have to use 
the binary complement of the random kink. 

We now calculate once more whether the first bond points at a negative or positive 
direction, and with this information we can calculate the new positions of the first monomers: 

''^i = '^i ~ 2«before + ^"after- \°) 

When the last monomer is selected, it can simply be replaced with random new bonds: 

yL-2 = y (9) 

ZL-2 = Z. 

E. Generation of random kinks 

The algorithm described above relies on the availability of random kinks. These kinks 
should be generated with the probabilities as given in equation (^. Since the two bonds in 
a kink have opposite directions only one bond has to be generated; the bond on the other 
side of the monomer is easily derived. 

The properties of detailed balance are used to create those bonds correctly. Certain 
properties must be enforced: first of all the x, y and z directions should occur with the same 
probability; secondly the ratio of the probabilities for + and - bonds is given by quotient of 
P^ and P+, as given in equation (^); this quotient is given by 

The first property is enforced by rotating some of the bonds (we used 50%) the following 
way: xn^y, yi-^z and zi— i>x. Using the randomly generated bit pattern r the following 
statements are used to rotate the bonds: 



X = (r A x) V (-ir A y) 

y = {r Ay)y (-ir A z) (11) 

z = {r A z) y (-ir A x). 

The second property is then enforced by inverting some of the bonds. With 50% prob- 
abihty, the negative bonds are inverted, with P^^^ times 50%, also the positive bonds are 
inverted. To make sure that all random kinks are independent we create a list of those and 
reshuffle this list regularly. 

III. SIMULATIONS 



The simulation algorithm described in sections |IIQ to ^l^ was implemented using the C 
programming language. The random number generator we used is a lagged (24, 55) additive 
Fibonacci generator. The simulations are done on a Silicon Graphics Origin 200 (180 Mhz) 
and on a DEC Alpha (466 Mhz) computer. The latter is faster and takes about 1.1 /xs per 
Monte Carlo step. 

The polymers where initialized in a U-shape with both ends in the direction of the electric 
field. At regular intervals we checked whether a polymer has moved at least its own size, 
which is the maximum distance between any two monomers. When this has occurred for a 
polymer, we assume that the polymer has thermalized, the real measurement starts when 
this thermalization has finished. 

The measurement is stopped when all polymers have thermalized and the average dis- 
tance traveled by all polymers is a few times their own size. We assume that measurements 
are statistically independent when a polymer has traveled a distance equal to its own size. 

We have performed simulations for lengths 3 up to 200. The time taken to calculate the 
drift velocity varied from a few seconds for small polymers up to about 17 hours for the 
longest polymers {L = 200) in the smallest electric field {E = 0.001). Simulations of longer 
polymers take too much time to compute the drift velocity for small electric fields. 



IV. SCALING ARGUMENTS FOR THE DRIFT VELOCITY 

The velocity of a polymer in a small electric field behaves according to the Nernst- 
Einstein relation, v = FD, where F = qLE is the force. The diffusion constant can thus 
be calculated from the drift velocity hj D = v/{qLE) in the limit E ^ 0. De Gennes ||l| 
found the diffusion constant to be proportional to D ~ L^^. This means the drift velocity 
is: V ~ qE/L. 

For slightly larger electric fields the Nernst-Einstein relation breaks down. Barkema, 
Marko and Widom [[7[ give an intuitive explanation of the dependence of the drift velocity 
on the electric field strength. The argument goes as follows. 

A random polymer will have an end-to-end length around h = yL. When an electric field 
is applied, the polymer is stretched in the direction of the electric field. When the electric 
field exceeds a certain level, the polymer as a whole does no longer resemble a random walk: 
h > ^/L. One may cut the polymer into Uf, pieces ('blobs') of length Lf, = L/rif,, that 
each still look like a random walk; the average end-to-end distance of the blobs is equal to 
(hf,) = v^II. 

Two forces work on the blob. The elastic force tries to contract a stretched polymer and 
is proportional to the size of the blob and inversely proportional to the length of the part 
of the polymer that forms the blob: -Feiastic ~ hh/Lf,. The electric force tries to stretch the 
polymer and is proportional to the size of the blob as well as the electric field: -Feiectric ~ hbE. 
These two forces have to be in balance which implies that the blob size is Lb ~ E~^. 

The Nernst-Einstein relation now applies to the blobs, so f = FbDb = qLbEDb- Again, if 
the blob size is large enough, D^ ~ L^^ which makes the speed of the polymer quadratic in 
the electric field: v ~ E/Lb ~ E'^. This effect has already been observed in the Rubinstein- 
Duke model by Barkema, Marko and Widom 0. 



V. RESULTS 

The results of our simulations are presented in figure |^. The short polymers, up to 
length 20, show a behavior different from the longer polymers. These short polymers have 
no superlinear dependence on the velocity on the electric field. 

When a small force, EL <^ 1, is applied to the polymers, the velocity of the polymers 
varies linearly with the electric field. When a force around EL ^ 1 is applied to the longer 
polymers, the polymer velocity depends superlinearly on the electric field. We show in figure 
^ that the dependence becomes quadratic for long polymer chains, as derived in section |T|. 
For much larger electric fields, the velocity decreases to zero. 

For polymer length L = 100 we performed some short simulations to get insight in the 
typical movement of the center of mass of the polymer. In figure ^ the position of the center 
of mass, scaled with a factor of l/E, is plotted as a function of time, for different field 
strengths. The starting point of the polymers are chosen in a way that the graphs do not 
overlap. For the smallest electric fields the movement is just like one would expect from a 
diffusing particle, it moves randomly, but with some preferential direction. For the electric 
field in the middle range, the diffusion effect becomes relatively smaller. This results in a 
smoother behaviour. In high electric fields the movement of the center of mass sometimes 
halts, when the force on the ends of the polymer pulls the polymer into a U-shape. When 
this happens the polymer has to untangle itself before its center of mass can move forward 
again. 

A. Polymer shapes 

The polymer shape in a small electric field resembles a random walk, as shown on the left 
side of figure ||. When the electric field is increased, the shape becomes stretched parallel 
to the electric field [|ll|]; the configuration may be viewed as a set of blobs which move with 



independent speed, as discussed in section IV. As shorter polymers move more quickly in 
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a given electric field, the blob- configuration moves faster than a random walk configuration 
which results in a superlinear increase of speed when the electric field is changed. 

When the electric field is increased above a certain value the shape may transforms into 
an U-shape, as shown in figure ^. With higher electric fields it becomes more difficult to 
escape from this U-shape. Since the polymer cannot move sideways it is trapped in the 
lattice for a long time compared to the time it moves. 

Figure ^ shows polymers in different configurations. The first polymer is stretched in the 
direction of the electric field. This configuration may be viewed as a large number of very 
small blobs. As such, the polymer has a high velocity, which may also be seen in figure § 
near 5.8 ■ 10'^ Monte Carlo steps. 

The second polymer is a transition configuration between the fast-moving cigar-like con- 
figuration as described above and the U-shape configuration. The polymer forms a 'knot' 
which is later passed by the trailing end of the polymer. 

The third polymer has a typical U-shape. The polymer in this configuration is almost 
fully stretched. This means that it has only a small number of kinks, which means that 
there is almost no stored length. 

Just before the polymer escapes from the U-shape, as is the fourth polymer, its con- 
figuration is very much stretched and has almost no stored length. This state transforms 
quickly in a state that resembles the state of the first polymer in figure |^. 

As described before, the number of kinks in a polymer decreases when the polymer is 
stretched. To check the dependence on the electric field we have performed some short 
simulations to find the average number of kinks on each location along the polymer. The 
simulations consisted of 10^ Monte Carlo steps after 2 ■ 10^ steps of thermalization, starting 
with a random configuration. Every 10^ Monte Carlo steps the kinks are counted. The 
fraction of time that a kink exists on a certain location is displayed in figure |^. 

For small electric fields the polymer configuration is known to resemble a random walk 
in three dimensions. The average number of kinks is thus expected to be 1/6. For higher 
electric fields the U-shape configuration becomes more frequent. In this configuration the 

11 



kinks are likely to diffuse towards the ends of the polymer, which means that the average 
number of kinks in the middle of the polymer decreases. When this happens we can no 
longer apply the blob argument as described in section |IV|. The mobility of the blobs in the 
middle of the polymer decreases as the average number of kinks in that region decreases. 

For a longer U-shaped polymer it becomes more difficult to escape from this configu- 
ration. The probability of a kink moving from one end of the polymer to the other end 
decreases with the length of the polymer. This means that longer polymers spend a longer 
time in the U-shape configurations. 

When the density of kinks becomes less than 1/6 per monomer, the entropic force that 
contracts the polymer is no longer in balance with the electric force. The polymer itself 
now transports the force along the chain, which may be better explained by the continuous 



model of Deutsch and Madden 12 



B. Comparisons to previous reports 

The results of the Duke-Rubinstein model have been compared to actual experiments 
13| . For longer polymers, the data is well described by 

1/2 



L'v 



a 



'leV (leV 



(12) 



This function is equivalent to the function v"^ = aE"^ + bE^, where a and b are functions of a, 
j3 and L. To check whether our results show the same scaling behaviour, we collapsed our 



data to the function v' = y E'"^ + E'"^ in figure ||, where v' = {Vb/a)v and E' = {Jb/a)E. 
The data in the third region is discarded in calculating a and b. 



Experiments have been performed on DC electrophoresis [0,|I^, both articles confirm 
the existence of regions where v ^ E and v ^ E'^. 

The diffusion constant is equal to D = ^JajL. The scaling found by Barkema and 
Krenzlin |jlO[| is given by DN^ = 0.173 + 1.9A^^^/'^, where N = L — 1. Figure ^ shows our 



results compared to their scaling function. Our results agree within statistical errors. 

12 



VI. CONCLUSIONS 

The cage model is extended to simulate DC electrophoresis, and the drift velocity of 
polymers in a gel is measured as a function of polymer length L and electric field strength 
E. 

The polymers behave differently in three regimes of the electric field: in a small electric 
field the velocity depends linearly on the electric field, in a high electric field the polymers are 
likely to be trapped in an U-shape. The regime in between shows a superlinear dependency 
on the electric field. We showed that this dependency becomes quadratic in the electric field 
for L —^ oo. 

We have shown that the average number of kinks is not equally distributed over the 
polymer in high electric fields, because the polymers tend to get stretched in the direction 
of the electric field. Hence the stored length disappears from the polymer. 
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FIGURES 
FIG. 1. Movement of a defect along the chain. When a defect moves along the chain, it displaces 

monomers which it passes by a distance equal to the stored length. 




FIG. 2. One elementary move of a monomer: a 'kink' (pair of anti-parallel neighbouring bonds) 
is replaced by another one. 





FIG. 3. Velocity of polymers of lengths 3 up to 200 in electric fields between 0.001 and 1. The 
velocity is plotted on the vertical logarithmic axis. On the horizontal logarithmic axis the electric 
field E is plotted. 
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FIG. 4. The position of the center of mass of a polymer of length 100 divided by the electric 
field for the electric fields 0.1, 0.03, 0.01, and 0.003. The lines are the expected positions using the 
average speed measured by the long experiments. 
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FIG. 5. Three polymers of length 100 in different electric fields. From left to right: E = 0.003, 
E = 0.01, E = 0.03. Polymers in small electric fields look like random walks. In slightly larger 
electric fields the ends tend to protrude. In high electric fields the polymer does not look like a 
random walk. 
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FIG. 6. Four polymers of length 100, and E = 0.1. From top to bottom the times in Monte 
Carlo steps are: 5.8 • 10^ 8.6 • 10^ 1.25 • 10® and 1.66 • 10^ 
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FIG. 7. The average number of kinks on monomer number. The polymers are of length 100 
and electric fields are 0.003, 0.01, 0.03 and 0.1. The line gives the expected value 1/6 of kinks in a 
random walk. 



0.2 



0.15 




0.05 



10 20 30 40 50 60 70 80 90 100 
monomer position 



FIG. 8. Transition between the linear and quadratic dependence of the velocity on the electric 
field. For various polymer lengths, the scaled velocity v' = {vb/a)v is plotted as a function of 
scaled electric field E' = {yJh/a)E, where a and b are L-dependent parameters given in table ||. 



The curve is given by v' = y E'"^ + E'^. 
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FIG. 9. Diffusion constant calculated from our measurements, compared to the scaling relation 
found by Barkema and Krenzlin. This scaling relation is a straight line when N'^D is plotted as a 
function of N^"^/^. 
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TABLES 
TABLE L Encoding of a bond in three bits, where x*-*' is the z*^ bit of x and so on. Note that 
the encoding of the negative bonds is the binary complement of the positive bonds. 



bond 


^{i) 


y« 


z(») 


+x 


1 








+y 





1 





+z 








1 


-X 





1 


1 


-y 


1 





1 


-z 


1 


1 






TABLE IL Values for a and 6, obtained by fitting the drift velocity to the form v"^ = aE"^ + bE ; 
these values are used for scaling in figure ^. For large L the parameter a decreases quadratically 
with L and b is more or less constant. The graphs of L^a and b show evidence of convergence to a 
constant; this is in agreement with Equation (^). 



L 


a 


b 


30 


1.85(6) • 10-^ 


1.9(2) • 10-2 


50 


4.4(2) • 10-5 


5.5(2) • 10-2 


70 


1.79(5) • 10-5 


8.3(2) • 10-2 


100 


7.32(7) • 10-*^ 


1.12(1) -10-^ 


140 


3.31(3) • 10-6 


1.29(1) • 10-1 


200 


1.5(1) • 10-*^ 


1.4(7) • 10-1 
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